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ABSTRACT 

The coalescence of a massive black hole (MBH) binary leads to the gravitational-wave recoil 
of the system and its ejection from the galaxy core. We have carried out iV-body simulations of 
the motion of a Mbh = 3.7 x 10 6 M MBH remnant in the "Via Lactea I" simulation, a Milky 
Way sized dark matter halo. The black hole receives a recoil velocity of Vkick = 80, 120, 200, 
300, and 400 kms -1 at redshift 1.5, and its orbit is followed for over 1 Gyr within a "live" host 
halo, subject only to gravity and dynamical friction against the dark matter background. We 
show that, owing to asphericities in the dark matter potential, the orbit of the MBH is highly 
non-radial, resulting in a significantly increased decay timescale compared to a spherical halo. 
The simulations are used to construct a semi-analytic model of the motion of the MBH in a time- 
varying triaxial Navarro-Frenk-White dark matter halo plus a spherical stellar bulge, where the 
dynamical friction force is calculated directly from the velocity dispersion tensor. Such a model 
should offer a realistic picture of the dynamics of kicked MBHs in situations where gas drag, 
friction by disk stars, and the flattening of the central cusp by the returning black hole are all 
negligible effects. We find that MBHs ejected with initial recoil velocities Vkkk ft 500 kms -1 do 
not return to the host center within a Hubble time. In a Milky Way-sized galaxy, a recoiling hole 
carrying a gaseous disk of initial mass ~ A/bh may shine as a quasar for a substantial fraction of 
its "wandering" phase. The long decay timescales of kicked MBHs predicted by this study may 
thus be favorable to the detection of off-nuclear quasar activity. 
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Introduction 



Intermediate-mass black holes may have formed at redshift z > 15 at the bottom of shallow dark- 
matter potential wells (jMadau fc Reed 120011 ). These seed holes may have grown through gas accretion 
and binary coalesce nces to become the supermassive variety that is ubiquitously found today at the center 
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of nearby galaxies (jKormendy et al 
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cold dark matter (CDM) cosmologies, where large halos are assembled via the hierarchical assembly and 
accretion of smaller progenitors, close MBH binaries inevitably form in large numbers during cosmic history 



(jBegelman et al 



1980; 



Volonteri et al.ll2003l ). The presence of a M BH binary with separation < 1 kpc ha s 



been revealed by Chandra observations of the nucleus of NGC 6240 (jKomossa et alJl2003t iMax et al.l l2007). 
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The Very Long Baseline Array (VLBA) discovery in the radio galaxy 0402+379 of a MBH binary 
system with a p r ojecte d separation of just 7 pc and a combined mass of ~ 1.5 x 10 s M Q was reported by 
Rodriguez et al. ( 20061 ). A MBH binary may shrink owing to stellar and/or gas dynamical processes (e.g., 



Maver et al 



20071 ) and finally coalesce when gravitational wave radiation dominates orbital energy losses. 



Recent developments in numerical relativity ( Pretoriusj|2005i Campanelli et al]|2006 ; Baker et al.|l2006a ) 



Herrmann et al. 


2007: 


Gonzalez et al. 


2007) 



Gonzalez et al.ll2007l ). Gravitational wave emission is typically anisotropic because of 
asymmetries associated with the masses and spins of the black holes, and causes the center of mass of the 
system to recoil in order to balance the linear momen tum carried away by gravitational radiation (jBekenstein 



1973c iFitchett fc Detweilerlll984l : iFavata et al.l 120041 ). The recoil velocity T4ick depends on the binary mass 
ratio <7b = M1/M2 < 1 on the dimensionless spin vectors of the pair Si and Hi (0 < en < 1), and on the 



orbital parameters. All current numerical data on kicks can be fitted by ( Baker et al. 2008h 
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where /j, = <7&/(l + qb) 2 is the symmetric mass ratio, 9i is the angle between the dimensionless spin vector 
Hi = Si/Mf of the ith black hole and orbital angular momentum vector, <pi is a projection angle between 
the spin vectors and a reference angle that lies in the orbital plane, and &i(qb) = $2(l/<7b) are constant for 
a given value of q b . Here, A = 1.35 x H^kms -1 , B = -1.48, H = 7540 ± 160 km s" 1 , £ = 215° ± 5°, and 
K = 2.4 ± 0.4 x 10 5 kms _1 . Assuming random spin orientations, qb > 1/4, and ai = 0,2 = 0.9, recoiling 



black holes can get a kick velocity > 500 kins 1 approximately 60% of the time (see Table 3 of Baker et al 
2008). For > 0.1, the percentage of kicks with > 500 kms -1 decreases to ~ 20%. Spins that are aligned 



with the orbital angular momentu m vector (as expected u nder the action of external torques provided by 
a circumbinary accretion flow, see iBogdanovic et al.l 120071 ) yield recoil velocities below 200 kms -1 , while 
the configuration producing the maximum recoil kick corresponds to equal-mass maxim ally rotating holes 
with a nti-aligned spins oriented parallel to the orbital plane, Vkick = Kfi 3 = 3750 kms -1 (| Campanelli et al 
2007al ). 



If not ejected from the host altogether, the recoil ing MBH will travel som e maximum distance and 
then return to the center subject to dynamical friction (jMadau fc Quataertil2004f) . Galaxy mergers are also 
a leading mechanism for supplying gas to their nuclear MBHs, and a recoiling hole can retain the inner 
parts of its accretion disk, providing fuel for a continuing luminous phase along its trajectory. Two possible 
observational manifestations of gravitat ional-radiation ejection have then been suggested: (1) spatially offset 
activ e galactic nuclei (AGN) activity (jMadau fc Quataertl |2004 iBlecha fc Loebl l2008t IVolonteri fc Madau 
2008); and (2) broad emission lines that are substa ntially shifted in velocity relative to the narrow- line 
gas left behind (jBonning. Shields, fc Salvianderl 120071 ). T he effect of gravitational wave recoil in the mass 
buildup of MBHs is more prominent at high redshifts (e.g., Volonteri fc Rees 20061 : Tanaka fc Haimanl 2009 ) 
and therefore the detection of offset nuclei is difficult. Observational evidence of re coiling MBHs is scarce 
and highly controversial. A recoiling SMBH candidate at z = 0.71 was reported by Komossa et al. ( 2008t ) 
in quasar SDSS J092712. 65+294344.0. The broad-line region of the quasi-stellar object (QSO), powered 
by a 6 x 10 8 A/ Q black hole, appeared to have a velocity offset of 2650 kms -1 with respect to the narrow- 
line region associated with the gal axy. However, sev e ral authors have challen ged this hypothesis, proposing 
that the object is a MBH binary (jDotti et alJl2008t IBogdanovic et al.l 120091 ) or an interacting galaxy pair 
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([Shields et al.ll2009l ; Eeckman et al.ll2009h 



The observability of recoiling MBHs depends sensitively on their dynamics in galaxy halos. The radial or- 
bit of a recoiling hole in a sp herically symmetric po t ential was first studied analytically bv lMadau fc Quataert 
([20041 ) and numerically by iBovlan-Kolchin et al.l (|2004 ). These early studies showed that large kicks (~ 
400 kms -1 ) can displace MBHs tens of kiloparsecs away from the center of a Milky Way-sized stellar bulge 
and that, after the kick, the MBH undergoes several oscillations before decaying back to the bottom of the 
potential. Most of the orbital energy is lost during the MBH passages through the center, where dynamical 
friction is most efficient: the cuspy central stellar density profile is flatte ned by the heating effect of dy nam- 
ical friction, and the MBH decay timescale correspondingly lengthened. iGualandris fc Merrittl (|2008l ) have 
recently substantiated these results by performing direct summation iV-body simulations of MBH recoil in 
spherical galaxies with binary-depleted cores. They found that initially the MB H loses its energy due to 
dynamical friction as predicted by Chandrasekhar's theory ( Chandrasekhar 1943h . When the amplitude of 
the motion has fallen to roughly the core radius, the MBH and core experience damped oscillations about 
their c ommon center of m ass, which decay until the hole reaches thermal equilibrium with the surrounding 
stars. IVicari et al.l ([20071 ) evaluated the effect of non-spherical galaxy geometries on kicked MBHs using 



triaxial models, and found significantly longer decay timescales than in equivalent spherical systems, as in 
a non-spherical potential the hole does not return directly through the dense center where the dynamical 
friction force is highest. Blecha & Loeb (2008) studied the trajectories of kicked holes in a two-component 
galaxy model that includes a spherical stellar bulge and a gaseous disk, and found that kicks with initial 
velocity Vkick 200 kins -1 in the plane of the disk are quickly damped out in t < 10 6 5 yr. 



In this paper, we revisit the problem using a different approach. We carry out full iV-body simulations of 
a recoiling MBH that is subject only to gravity and dynamical friction against the dark matter background, in 
a high-resolution, non-axisymmetric, "live" potential. The host is the main halo of the Via Lactea I (hereafter 
VL-I) cosmological simulation (Diemand et al. 2007a, 2007b). We follow the MBH orbital behavior starting 
at redshift z = 1.5 (when the kick is assumed to occur) for more than 1 Gyr, as the host grows in mass 
and changes its shape from prolate to triaxial. We show that, owing to departures from axisymmetry in the 
dark matter potential, the orbit of the hole is highly non-radial, resulting in a significantly increased decay 
timescale compared to a spherical halo. The simulations are used to construct a more realistic semi-analytic 
model of the motion of the MBH in a time-varying triaxial Navarro-Frenk- White (NFW) halo plus a fixed 
isothermal stellar bulge, where the dynamical friction force is calculated directly from the velocity dispersion 
tensor. Such a model should offer a more realistic picture of the dynamics of kicked MBHs in situations 
where gas drag, friction by disk stars, and the heating effect of the returning hole on the central cusp are all 
negligible. 



2. Simulations setup and properties of the host 



The VL-I simulation was performed with PKDGRAV (Stadel 2001) a cosmological tree code that in- 
cludes gravitational multipoles up to hexa-decapole order to reach high accuracy in the force c alculation. It 
empl oyed multiple mass particle grid initial conditions ge nerated with the GRAFIC2 package ([Bertschinger 
200 1 ) in a WMAP 3-year cosmology (jSpergel et. alll2007t ). A bug in the original GRAFIC2 code caused the 



power spectrum used for the VL-I refinements to be that of the baryonic component, equivalent to an effective 
spectral index of n = 0.90 instead of the intended 0.95. In this cosmology subhalo concentrations and peak 
circular velocities are slightly lower than in WMAP 3-year, while erg and the main halo properties remain the 
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samejj The high-resolution region was sampled with 234 million particles of mass m p = 2.1 x 1O 4 M and 
evolved with a force resolution of e — 90 pc. It was embedded within a periodic box of comoving size L = 90 
Mpc, which was sampled at lower resolution to account for the large-scale tidal forces. The host halo mass 
at z = is M200 = 1-8 x 10 12 M within a radius of -R200 = 389 kpc (defined as the radius within which the 
enclosed average density is 200 times the mean matter value). In this work we rerun VL-I using PKDGRAV 
from redshift z% = 1.54 to Zf = 1.15, and follow the orbits of all dark matter particles as well as a new 
MBH particle placed at the center of the host. As in the original VL-I simulation, we employ a gravitational 
softening of 90 pc for the dark matter particles and the MBH, as well as adaptive time steps as short as 
r = 68, 500 yr, sufficient to ensure convergence in the density profile down to a radius of r conv ~ 1.0 kpc and 
to accurately sample the orbit of the MBH. The time-stepping criterion is given by At < 0.2^/e/ai, where ai 
i s the local acceler ation. The resolution of VL-I allows us to adopt the mass of SgrA*, Mbh = 3.7 x 10 6 M Q 
(|Ghez et al.ll2005h . 'or the central MBH particle: this implies a MBH-to-particle mass ratio of 175, enough 



to accurately reproduce the effect of dynamical friction. 

Large kicks can displace MBHs sufficiently far away that their decay times become a significant fraction 
of the age of the universe. It is interesting to look at the evolution of the host halo in terms of it s time- varying 



spherically averaged density profile and shape parameters. The fitting formula proposed bv iNavarro et al 



(|1997I ) provides a reasonable approximation to the density profile, 



x(l + x) z 

where x — r/R s and R s is the scale radius. The mass profile is given by M(< x) = M 2 oo f (%) / f (c) , where 
f(x) = ln(l + x) — x/(l + x) and c = -R200/ Rs is the concentration parameter. The escape speed from the 
halo center is 

where V^qq = GA^oa/Raoa- The quantities p s , R s , R200, M200, Vmax (the maximum circular velocity of the 
host) and f CS c(0) are given in Table 1 at different scale factors, starting with the time when the kick is 
imparted. 

CDM halos are known t o show significant de partures from sphericity (for a recent summary, see Allgood 



et al. 2006). As detailed in iKuhlen et alJ |2007j), we approximate the shape of the VL-I host potential by 



diagonalizing the unweighted kinetic energy tensor 

Kij = l}^2 v i,nVj,n, (7) 
n 

where Kij is related to the potential energy tensor Wij — ^XidQ / dxj through the tensor virial theorem 

1*± = 1K II + W U . (8) 

Here, 

Iij= J2^IL. (9) 



1 Note that thi s problem does not affect the more recent "Via Lactea II" and "GHALO" simulations (Dic mand et alj 2008; 
IStadel et ahlboosl) . 
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and r n = \/( x2 i + (jj n /q) 2 + i z n/s) 2 )- We assume cPlij/dt 2 = so that the eigenvectors of K{j reflect the 
principal axes of the potential ellipsoid. The latter is signific antly rounder than the mass distribution, and 
neither its shape nor orientation vari es much with the radius (jKuhlen et al.ll2007l ). The degree of triaxiality 
of the halo potential, T, is given by (jFranx et al.lll991l) 



T 



1 



(10) 



where q — b/a and s = c/a are the time-dependent intermediate-to-major and minor-to-major axis ratios, 
respectively (a > b > c). A halo is said to be oblate for T < 1/3, triaxial for 1/3 < T < 2/3, and prolate for 
T > 2/3. Figure [1] shows the evolution of the potential shape parameters with redshift at different radii. In 
the inner regions the axis ratios remain approximately constant after around z = 0.8, but before z = 1 there 
are significant changes in the outer regions, as the halo becomes more spherical. The triaxiality parameter 
remains mostly in the prolate regime (> 2/3) in the inner regions, while in the outer halo evolves from 
prolate at z > 1 to triaxial or slightly oblate at 0.7 < z < 1, to back to prolate at later times. Note that the 
VL-I host accretes some fairly massive subhalos between z = 1 and z = 0.5. Dynamical friction causes these 
subhalos to spiral in to the center over a few orbits, and they lose most of their mass in this process. The 
associated redistribution of material probably contributes to the observed shape adjustments. 



3. Dynamics of recoiling holes 

3.1. Orbits in numerical simulations 

We placed the MBH particle at the position of the densest point of the main VL-I halo at an initial 
redshift Z{ — 1.54, 300 Myr after the last major merger. At this epoch the host has M200 = 1-02 x 10 12 M 
and -R200 = 187 kpc. The kick was oriented at an angle of 20° to the minor axis of the host halo at Zj. The 
MBH orbit was tracked at every time step in our simulations, and its position and velocity were measured 
with respect to the central position and center of mass velocity, respectively. The five halo+MBH runs 
—corresponding to kick velocities Vkick = 80, 120, 200, 300, and 400 kms" 1 and labeled VL080 to VL400— 
were evolved for 1.15 Gyr (i.e., until a final redshift Zf — 1.15). All kick velocities are below the escape speed 
at Zi, v csc (r = Q,Zi) = 488 kms -1 . Each run consumed 13,000 CPU hours on the Pleiades Supercomputer 
Cluster at UCSC, and followed the MBH for 10,000 time steps. 

The resulting trajectories are shown in Figure [2] the orbits' parameters are listed in Table 2, and the 
three-dimensional rendering of the orbits in simulations VL120, VL200, and VL300 are shown in Figure [31 
Only Vkick > 300 kms -1 trajectories actually sample the outer halo with pericenter distances i? max ^ 30 
kpc, and only Vkick < 120 kms -1 trajectories return within 0.5 kpc from the center during the duration of 
the simulation. The motion of the hole remains nearly rectilinear for one or two oscillations only, as the 
y- and z-components of its orbit become rapidly important due to asphericities in the halo potential. This 
increases the MBH decay timescale compared to a spherical model, as we show below. Dynamical friction 
has only a weak effect on the maximum displacement of the MBH. This can be seen in Figure [3] (right top 
panel) , where a sixth simulation was carried out with the recoiling hole treated as a massless test particle of 
initial kick velocity Vkick = 80 kms" 1 . A comparison with VL080 shows how, for the first 2-3 oscillations, 
dynamical friction does not strongly influence the motion of the hole, and the maximum displacement is 
similar to that of the energy-conserving orbit. It is only at later times that the effect of friction sets in, 
reducing the amplitude and period of the oscillations and bringing the hole back to the center. Note how, 
for Vkick > 120 kms -1 , and because of the aspheric nature of the halo, the MBH spends most of its time 
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> 0.8 kpc away from the center and does not have a significant dynamical heating effect on the dark matter 
distribution in the nucleus. 



3.2. Orbits in a spherical NFW halo 



It is interesting at this stage to compare the results of our numerical simulations with a semi-analytic 
model of the motion of a recoiling MBH in an NFW halo. Such a model will allow us to follow the trajectory 
of a recoiling black hole for a Hubble time or until it returns to the center. We define the return time, i re turn, 
as the time it takes for the MBH to decay to within r = 1 pc of the center of the halo with \E/E- m \ < 0.001, 
where E is the total energy (kinetic + potential) of the MBH and E- m is its initial energy. The energy 
condition is set to ensure that the MBH is not simply going through a close periastron passage. 

We start by approximating the potential as spherically symmetric and static, with the z — host halo 
parameters given in Table 1. Under these assumptions the trajectory is purely radial, and the damping force 
from the background dark matter can be approximated by the classical Chandrasekhar dynamical friction 
formula (|Chandrasekharlll943t iBinnev fc Tremaindll987l ). The corresponding equation of motion is 



dv 



GM{< r) _ AttG 2 In Ap(r) Af BH 



erf(X) 



2X 



-X' 



v, 



(11) 



where X = v/\/2a. 



The proper definition of the Coulomb logarithm, In A, has been extensively debated. It is generally 
defined as In (6 max /6 m i n ), where the maximum impact parameter 6 max is the scale radius R s of the dark 
matter distribution, and the mini mum impact parameter 6 m , n is the radius of influence of the MBH, -Rbh = 
GM B n/cr 2 . Several studies (e.g.. IColpi et alJll999t lHashimoto et alj 12003) have shown that a dynamically 
varying value for In A provides a better estimate of dyn amical frictio n than a constant value when compared 
to A^-body simulations. Here we follow the treatment of lMaoa (|l993l ). and in the approximation of a spherical 
NFW host write the Coulomb logarithm as 



In A 



por 



dx 



PO Jxo x 2 {l + x) 2 



2 In 



1 



2x + 1 



x(x + 1) 



(12) 



where po is the central mass density, x = r/R s , d can be interpreted as the minimum impact parameter of 
the Chandrasekhar formula. Throughout this paper we use po = p(r = 20 pc) and set d — Rbh, the radius 
of influence of the MBH. 

The one-dimensional velocity dispersion a for an NFW profile can be solved numeric ally from the Jeans 
equat ion or approximated analytically for x — r/R s between 0.01 and 100 by the function (jZentner fc Bullock 
20031) 



x(l + x 2 ) 



x' 3 (l + x')'- 



■dx' 



1.4393a: 



0.354 



1 + 1.1756a; - 725 



(13) 
(14) 



We integrate the equation of motion numerically using an adaptive Adams-Bashforth-Moulton integration 
scheme. The resulting radial orbits for kick velocities Vkick = 80, 120, and 200 kms -1 are shown in Figured] 
(left panels). The decay timescale of a recoiling hole in a spherical potential is significantly shorter compared 
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to the results of TV-body simulations, mainly because of the efficiency of dynamical friction during each 
passage through the nuclear regions. In the Vkick = 120 kms -1 case, for example, the MBH is back to the 
center after 0.6 Gyr in the spherical case, while it is still wandering close to i? ma x in the simulations. A 
self-consistent estimate of the decay timescale must include the flattening of the cuspy central density profile 
by the oscillating hole. Such a cumulative heating effect, however, is negligible in this case, since due to 
the triaxiality of the potential the MBH does not affect the central density and velocity dispersion profiles 
dramatically. 



3.3. Orbits in a triaxial NFW halo 



The next-order approximation is to model the motion of the recoiling hole in a triaxial, dynamically 
evolving NFW dark matter halo, using the VL-I halo parameters given in Table 1 and the potential shape 
parameters plotted in Figure [TJ The orbit of the hole is fully specified by the conservative force of the dark 
matter potential V<I> and the damping frictional term: 



^ = -V$ + / DF , 



where 



and 



$ = 



GM 200 ln(l + r e /fl 8 ) 
" /(c) 



— I ™2 



1/2 



(15) 



(16) 



(17) 



is the ellipsoidal radius. Here q and s are the time- and radial-dependent axis ratios defined in Section 2, 
and x,y,z are Cartesian coordinates along the principal axis of the potential ellipsoid. Equation (jlip is no 
longer valid in a triaxial system, where t he velocity dispers ion is non-isotropic and the velocity distribution 
deviates from Maxwellian. We adopt the 



to triaxial systems (see also lVicari et al 



2007) 



Pesce et al.l (jl992T ) generalization of the dynamical friction formula 



DF 



-T a V a e a - T b V b i b - T c V c e c , 



(18) 



where Vj are the components of the black hole velocity along the principal axes ej = {e a , e b , e c } of the local 
velocity dispersion ellipsoid with a > b > c, and Ti are the dynamical friction coefficients. These arc given 

b y _ 

_ 2V2^G 2 p(r e )lnA(M BH ) 



x Bi(V, a), 



(19) 



where the velocity dispersion integral is given by 



Bi 



oo expC-XV^-^-) i 



■ du, 



(20) 



ti = <7j/eri, of — {o'a>' J b' a c} i s the velocity dispersion along the direction {e a ,eh,e c }, <j\ is the largest 
eigenvalue, and p(r e ) is the local mass density at the MBH's elliptical radius. In order to calculate the 
triaxial density profile, we deform the spherical density contours in such a way that the volume is preserved. 
In this approximation, the characteristic elliptical radius of the halo becomes i? e ,200 = (<7 s) -1 / 3 i?200- 
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A correct estimate of the velocity dispersion as a function of radius and redshift is crucial in the 
calculation of dynamical friction. Here, we take the following approach. First, we measure the "true" shape 
and orientation of the local velocity dispersion ellipsoids directly from the VL-I simulation (model A; see 
Section 13.3.11 for details). Next, we construct a model to calculate the velocity dispersion from the Jeans 
equation (model B; see Section f3. 3. 2| . In model B we neglect streaming motions and assume that the local 
velocity dispersion ellipsoids are aligned with the global potential shape, which results in an overestimate of 
the velocity dispersion integral Bi given by Equation [20] To normalize model B to the fiducial model A, we 
introduce a linear fitting factor r\ where Bf lodclA = rjBf lodcl B . The main characteristics of the MBH orbits 
are well reproduced by models A and B for a large range of recoil velocities using r\ = 0.5 (see Figure [5]). 

The resulting orbits are shown in the right panels of Figure HI The triaxial halo model qualitatively 
reproduces the results of the simulations, the highly non-radial MBH trajectories, and the extended wan- 
dering times of kicked holes. Return timescales exceed 10 Gyr already for Vkick = 200 krns -1 (see the last 
column of Table 2). 



3.3.1. Model A: Local Velocity Dispersions Measured in VL-I 



The local properties of the halo relevant to the calculation of dynamical friction were measured from the 
V L-I simulation as a function of redshift for 10 snapshots in the range < z < 1.54, following the method 



Zemp et al.l (|2009l ). At each redshift seven distances r = 1, 8, 25, 50, 100, 200, 400 kpc from the halo center 



of 

were randomly sampled with 10 spheres of radius 



r S ph(r) 



r sp h( 



;k P c) 



p(8kpc) 
p{r) 



1/3 



(21) 



where r sp h(8kpc) = 0.5 kpc and p(r) is the spherically averaged mass density at radius r. In each sphere 
we measure the local density and calculate the six components of the symmetric velocity dispersion tensor, 
o~ij = (viVj) — (vi) (vj) (here the indices i and j indicate the components along the principal axes of the global 
potential ellipsoid). We then diagonalize the dispersion tensor to obtain a set of eigenvalues and eigenvectors. 
The eigenvalues, {o^, a^, tr^}, are the components of the velocity dispersion in the ii basis. 

For computational convenience, we fit an analytical function to the mean valu e of the local velocity 
dispersion in all spheres at each radii. This function has the form (jPesce et al.lll992f) for model A: 



1 + Bir r ' 



1 + Dir 71 



-r/Ci 



(22) 



where Ai through Di and rrii,ni are the best fit values to the velocity dispersion profile in the ith direction 
at a given redshift. The parameters at z — are given in Table [3] and the corresponding best-fit curves for 
<Ji^A are shown in Figure [5Jd. 

The orientation of the local velocity ellipsoids with respect to the global shape was also measured as 
a function of radius and redshift. Table [4] shows the angles between the major, medium, and minor axes 
of the velocity dispersion ellipsoid and their counterparts in the global potential ellipsoid (a,/3, 7 respec- 
tively) averaged over the ensemble of spheres. The principal axes of the velocity ellipsoid show significant 
misalignment with the principal axes of the global potential shape: the distribution of orientation angles is 
quite isotropic and cannot be fit by a simple function. In our fiducial semi-analytical model (model A), the 
orientation of the local velocity dispersion is obtained by interpolating a grid of mean orientation angles as 
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a function of position and redshift at each time step of the numerical integration. Then a random value is 
drawn in the range allowed by the dispersion associated with the mean. 

The local density profile is shown in Figure [5^. The points represent the average local measurement, p, 
and the error bars are the dispersion around p, labeled a(p) in Table [4] The solid line represents the best 
fit NFW profile to the local density average at z = (see Table HJ. 



3.3.2. Model B: A simple treatment of local velocity dispersion 

While our fiducial model accurately reproduces important features of the orbits of MBHs in a tri- 
axial potential, having a simple prescription to calculate the velocity dispersion analytically would allow 
us to generalize our model and include the effect of other galactic components (see below). In this toy 
model, we assume that the local velocity dispersion ellipsoids are aligned with the potential shape: therefore 
{e Q ,eb,e c } = {e x ,e y ,e z } and all off-diagonal terms of the local velocity tensor vanish. We further assume 
that the halo is in steady state at each snapshot and that there are no streaming motions. Under these 
assumptions we solve for the velocity dispersion along the ith coordinate from a simplified Jeans equation: 

&i,B = ~~ r P(re)^A *4, (23) 

Pe J Xi ox i 

where p e is the density at the elliptical radius corresponding to the position of the MBH. We normalize the 
velocity dispersion integral (Equation to nBi(V, erg) in order to match the results of model A. Figure [5] 
shows a comparison of models A and B: maximum displacement distance and return times are accurately 
reproduced by model B for a large range of kick velocities with 77 = 0.5. This analytical representation of 
the velocity dispersion in a triaxial potential proves useful in the construction of the composite potential 
described below. 



3.4. Orbits in a triaxial NFW halo plus a stellar bulge 

A realistic study of the trajectories of recoiling holes must include the gravitational and frictional effect 
of a stellar bulge. Our final set of semi-analytic orbit integrations uses a two-component galaxy model 
consisting of a time-varying triaxial halo (with same parameters as above) and a fixed spherical bulge of 
stellar density 

P * ir) = 2nG(rhrty ™ 

with isotropic stellar velocity dispersion a* = 75 kms -1 , suitable for a Milky- Way-sized host. In the inner 
regions of the bulge, where stars are the dominant source of dynamical friction, the sphere of influence of 
the black hole is given by Rbh = GMbh/c»- The stars within this radius are bound to the black hole and 
do not exert dynamical friction, and therefore a MBH traveling through the very center of the bulge will 
experience an effective core radius r c = -Rbh- We truncate the bulge profile at an outer radius of r\> = 3 kpc 
in order to obtain a finite bulge mass at large radii, where the dark matter halo dominates the potential. In 
this model, the mass of the stellar bulge within the outer truncation radius is M*(< r^) = 8 x 10 9 M Q . 

To find the velocity dispersion tensor of the composite profile, erf.-, we solve the Jeans equations under 
the assumption that the velocity ellipsoid is aligned with the axes of the (dynamically evolving) triaxial 
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NFW potential. Thus, the three principal components of the velocity dispersion tensor are given by 

i r ao> tot 

CT i,tot = / Ptot -5 dXi (25) 

Ptot J Xi OXi 

where ptot and <&tot are the total (NFW halo + stellar bulge) density and potential. We calculate the Coulomb 
logarithm from Equation [T2l using the total mass density and d — GMbk/ct 2 , where the composite velocity 
dispersion is now given by er tot = \J a\ + a 2 + <r 2 (with Oi given by Equation |25|) . As in Section [3.3.21 we 
normalize the velocity dispersion integral to r]Bi(V, <?t t), where Bi(V, a tot) is the velocity dispersion integral 
of the composite potential. We assume that the spherical stellar bulge fully dominates the potential in the 
region r < 100 pc and therefore dynamical friction is well approximated by the Chandrasekhar formula. 
We fit for 77 by comparing orbits obtained with Equation [TT] with those obtained with Equation [18] for 
Vki C k < 250 kms" 1 . The best fit yields r\ = 0.1. 

Table 4 gives the MBH apocenter, its pericenter, and the return time calculated using our two-component 
model: (1) for Vkick > 460 kms" 1 we stopped numerical integration after a Hubble time <h, while the hole 
was still wandering tens to hundreds of kiloparsecs away from the center; (2) for kick velocities below 
380 kins" 1 , dynamical friction against bulge stars now efficiently damp the motion of the MBH already on 
the first outward trajectory, and reduces the decay timescale to less than 2 Gyr. Recoiling holes do not leave 
the bulge; (3) for the maximum kick velocities predicted in the case of non-rotating holes, Vkick ^ 200 kms -1 , 
the MBH reaches a maximum distance of only 40 pc from the center and decays back within 2 Myr; (4) black 
holes that leave the stellar bulge and enter the triaxial dark matter halo do not return to the center within 
a Hubble time. The pericenter distances, apocenter distances, and the return times of MBHs are shown in 
Figure [7] for a dark matter only potential and a more realistic dark matter + bulge potential. According to 
the latter model, a MBH which is kicked with initial velocity Vkick = 400 kms" 1 reaches i? max before 10 8 yr, 
a time comparable with the typical QSO lifetime, and spends most of its time orbiting at a distance r > 1 
kpc away from the center of the bulge. 



4. Summary 

Coalescing MBH pairs will give origin to the loudest gravitational wave events in the universe, and are 
one of the primary targets for the planned Laser Interferometer Space Antenna (LISA; e.g., Sesana et al. 
2004) . The anisotropic emission of gravitational waves also removes net linear momentum from the binary 
and imparts a kick to the center of mass of the system. The outcome of this "gravitational rocket" has been 
the subject of many recent numerical relativity studies. Non-spinning holes recoil with velocities below 200 
kms -1 that only depend on the binary mass ratio, while much larger kicks require rapidly rotating holes. 
Little is known about the masses of MBH binaries and their spins: the distribution of all binary mass ratios 
expe cted in some hierarchical models of the co-evolution of MBHs and their hosts is found to be relatively 



flat (jVolonteri fc Madaul 120081 ): if it is not ejected from the host altogether, the recoiling MBH will travel 



some maximum distance and then return towards the center on a decay timescale that depends on the shape 
of the potential and on the effectiveness of gas drag and dynamical friction against the stars and the dark 
matter of the host galaxy. 

We have carried out a detailed study of the fate of bound recoling holes in Milky Way-sized potentials, 
running TV-body simulations of the motion of a Mbh = 3.7 x 10 6 M Q MBH remnant in the "Via Lactea I" 
dark matter halo. In the simulations, the MBH receives a kick velocity of Vkick = 80, 120, 200, 300, and 400 
kms" 1 following the coalescence of its progenitor binary, and moves within the "live" host subject only to 
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gravity and dynamical friction against the dark matter background. We have used these calculations to build 
realistic semi-analytic models of the hole's trajectory in a time-varying triaxial NFW potential, where the 
dynamical friction force is calculated directly from the velocity dispersion tensor, and in a two-component 
triaxial halo+spherical bulge model. The latter case should offer a more realistic picture of the dynamics of 
kicked MBHs in situations where gas drag, friction by disk stars, and the heating effect of the returning hole 
on the central cusp are all negligible. Our results on the trajectories of recoiling MBHs can be summarized 
as follows: 



Owing to asphericities in the dark matter potential, the black hole's orbits are highly non-radial, result- 
ing in a significantly increased d ecay timescale comp ared to the spherical case. This is in qualitative 



agreement with earlier results bv lVicari et al.l (|2007l) 



2. In a triaxial NFW halo return timescales to the center exceed 5 Gyr already for Vkick = 200 kms 1 , 
and are longer than the Hubble time for Vkick > 250 kms -1 . 

3. In a triaxial halo+spherical bulge potential, decay timescales are much shorter than in the bulgeless 
case. For kick velocities Vkick < 380 kms -1 , dynamical friction against bulge stars now efficiently 
damp the motion of the MBH already on the first outward trajectory, and reduces the decay timescale 
to less than 2 Gyr. For recoil velocities Vkick > 500 kms -1 the MBH does not return to the center of 
its host within a Hubble time. Recoling black holes do not leave the bulge and remain within a few 
tens of parsecs from the center for Vkick 200 kms -1 . 



A kicked MBH can retain the inner parts of its accretion disk, providing fuel for a continuing luminous 
phase along its trajectory. Let us assume all recoiling holes accrete at a fraction Je of the Eddington rate 
Me = 4:TrGMBnrn p /(caT£), where e is the radiative efficiency. The duration of the luminous phase depends 
on the amount of disk material out to the radius f? ut ~ C^bh/V^^ that is carried by the hole. In the case 
of an a-disk, this is given by (Loeb 2007) 

Af disk « (1.9 x 10 6 Mq) a :;- 8 ( e _ 1 // £ )-°- 6 M 7 2 - 2 V 3 - 2 - 8 , (26) 

where e_i = e/0.1, M7 = A/bh/IO 7 Mq, V3 = Vki c k/10 3 kms -1 , and a_i = a/0.1 is the viscosity parameter. 
The condition Afdisk < A^bh then requires 

V^k > 550 kms -1 aI°- 2 Vi//Er - 21 ^7 ' 43 . (27) 

For lower kick velocities Afdisk = Mbb., corresponding to an AGN lifetime of £qso = eccry/ (4wGm p fE) ~ 
4.5 x 10 7 yr (e-i/Zs). A recoiling hole/disk system with (Mr, a, e, fE, Afdisk) = (1, 0.1, 0.1, 1, Mbh) could 
then be shining for half a Gigayear as an off-center quasar over a large fraction of its "wandering" phase. 
Thus, cases where the recoil kick is large enough to launch the MBH into the triaxial halo are favorable for 
the detection of off-nuclear quasars. However, if the MBH is initia lly embedded in a gas-rich environment, 
gas drag may damp its motion significantly ( Guedes et al. 2008h . even for moderate kicks, lowering the 



detection probability. Furthermore, the spins of both black holes in a MBH binary te nd to align due to 



torqu es induced by the surrounding gas, reducing the kick velocity to Wkick < 200 kms 1 (jBogdanovic et al 



2007). The motion of a recoiling MBH in a gas-rich merger including a stellar and dark matter component 



will be the subject of a subsequent paper. 
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Fig. 1. — Host halo (potential) shape parameters as a function of redshift at different ellipsoidal radii. Intermediate-to-major 
axis ratio q {solid points), minor-to-major axis ratio s {empty circles), and triaxiality parameter T {insets). 

Insets have the same x-axis range as the main plots. 
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Fig. 2. — Response of a M B h = 3.7 X 10 6 Mq MBH to a kick at Zi = 1.54 in the aspherical potential of the "live" VL-I Milky 
Way-sized halo. The radial distance R of the hole from the center is plotted vs. time for Vkick = 80, 120, 200, 300, and 400 
kms" 1 . Each orbit was sampled with 10,000 points. 



0.0 0.1 0.2 0.3 0.4 

Time [Gyr] 




Fig. 3. — Left top, left bottom, and right bottom panels: Three-dimensional orbits of the recoiling MBHs in the VL120, VL200, 
and VL300 simulations. The first 0.5 Gyr are plotted in yellow, the following 0.5 Gyr in red, and the remaining 0.1 Gyr in 
purple. Box sizes are 2.5, 7.6, and 25.1 kpc, respectively. Right top panel: Comparison between orbits in VL080 (yellow) and 
the corresponding energy-conserving orbits in the massless hole simulation (green). 
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Fig. 4. — Decay histories of recoling MBHs. The radial distance R from the center is plotted as a function 
of time for a spherical NFW (left panel) and a triaxial NFW halo (right panel). The ./V-body simulation 
results (red curves) are superposed to the semi-analytic orbit integrations according to model A (orange). 
The insets are a close-up of the respective orbit over a timescale of 1.1 Gyr. 
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Fig. 5. — Properties of the VL-I main halo at z — 0. Left: Local density averaged over an ensemble of spheres 
at discrete radii from the VL-I halo center (asterisks) and best fit NFW profile (solid line). Right: Average 
local velocity dispersion along the principal axes of the local velocity dispersion ellipsoids as a function of 
radius (asterisks) and best fit velocity dispersion profile (solid line). The error bars represent the dispersion 
around the mean value. 




Fig. 6. — Left: Maximum displacement distance in model A (fiducial) compared to model B for kick velocities 
in the range 80 < Vkick < 600 kins" 1 (asterisks). Right: Return times of models A and B for kick velocities 
80 < Vkick < 250 kms -1 (asterisks). Colors represent magnitude of the recoil velocity from Vkick = 80 kms -1 
(blue) to Vkick = 600 kms^ 1 (red). 
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Fig. 7. — Upper panel: A set of apocenter distances (solid line) and pericenter distances (dashed line) for a 
recoiling MBH of mass M. = 3 x 10 6 M in a triaxial Milky Way-sized dark matter only host (green) and 
dark matter + bulge host (blue). The colored areas show the corresponding regions in the R — Vkick plane 
occupied by the wandering holes. Lower panel: Return timescales of a MBH in a dark matter only host 
(green line) and a dark matter + bulge potential (solid blue line). Also shown is the time it takes to the hole 
to reach its apocenter (dashed- dotted blue line). 
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Table 1: VL-I halo parameters 



a 


Ps 

(10 6 M Q kpc" 3 ) 


R s 
(kpc) 


-R200 
(kpc) 


M 200 
(10 12 M ) 


^max 

(km s -1 ) 


^esc 

(km s _1 ) 


0.393 


0.16 


38.0 


194.3 


1.03 


160.53 


488.5 


0.423 


0.21 


36.7 


213.1 


1.13 


163.7 


498.2 


0.465 


0.41 


31.3 


233.8 


1.19 


167.9 


510.9 


0.507 


0.54 


30.8 


250.9 


1.22 


166.7 


507.3 


0.549 


0.72 


29.7 


271.5 


1.31 


170.4 


518.4 


0.591 


0.99 


28.1 


292.5 


1.43 


176.4 


536.7 


0.633 


1.40 


26.2 


311.8 


1.54 


182.6 


555.8 


0.675 


1.87 


25.1 


327.4 


1.61 


186.4 


567.2 


0.762 


2.40 


26.0 


356.6 


1.76 


189.2 


575.6 


0.877 


3.54 


26.2 


376.2 


1.77 


187.0 


569.0 


0.901 


3.67 


26.7 


379.2 


1.77 


185.6 


564.8 


0.950 


4.51 


26.2 


384.9 


1.77 


185.9 


565.5 


1.000 


5.33 


25.8 


389.3 


1.77 


185.1 


566.2 



Table 2: N-body Simulation Results 



Run Name 


Hick 

(km s _1 ) 


^max 

(kpc) 


(kpc) 


^-end 

(Gyr) 


^end 

(kpc) 


^return 

(Gyr) 


VL080 


80 


1.18 


0.03 


1.15 


0.09 


1.16 


VL120 


120 


2.49 


0.56 


1.15 


1.90 


2.78 


VL200 


200 


7.69 


0.72 


1.15 


3.71 


8.45 


VL300 


300 


28.21 


1.51 


1.15 


14.93 


> t u 


VL400 


400 


83.65 


22.95 


1.15 


22.95 


> t n 



Columns 2,3,4,5,6 and 7 list the initial kick 
velocity, the MBH apocenter, its pericenter, the 
end time of the simulation, the distance of the 
MBH from the halo center at t cn d, an d the 
return time calculated using a triaxial NFW 
model, respectively. The return time, i rc turn, is 
defined as the time is takes for the MBH to lose 
all but 0.1% of its initial total energy and decay 
to within 1 pc of the center of the halo. 
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Tabic 3: Best- fit parameters to the velocity dis- 
persion profile at z = 0. 





A 

(km 2 s- 2 ) 


B 

(kpc-" 1 ) 


C 
(kpc) 


D 

(kpc-™) 


m 


n 


°i 


2.24 x 10 4 


1.145 


172.21 


0.0026 


-4.87 x 10 3 


1.1132 


-I 


7.16 x 10 2 


0.567 


153.55 


14.300 


-1.29 x 10 2 


0.1217 


°l 


3.65 x 10 2 


0.4102 


117.02 


22.698 


-0.19 x 10 2 


0.1646 



Table 4: Summary of Local Properties at z = 



r 


(kpc) 


1 


8 


25 


50 


100 


200 


400 


P 


(M pc- 3 ) 


7.90xl0- 2 


6.65xl0- 3 


5.64xl0- 4 


2.80 xlO- 4 


3.35 xlO- 5 


3.65 xlO- 5 


1.76 xlO" 6 


&a 


(kms -1 ) 


102.1 


148.0 


143.7 


148.1 


126.9 


121.1 


79.31 


&b 


(kms- 1 ) 


83.38 


125.1 


124.5 


116.4 


85.98 


79.60 


50.06 


&c 


(kms- 1 ) 


77.99 


113.4 


116.4 


106.3 


77.26 


63.54 


38.71 


a 


(°) 


19.82 


53.19 


61.45 


62.30 


36.57 


34.79 


45.19 


B 


(°) 


69.17 


48.28 


57.47 


50.43 


54.81 


54.66 


57.69 


i 


(°) 


84.10 


69.04 


69.31 


66.32 


54.63 


58.78 


56.22 


<p) 


(M pc- 3 ) 


5.88xl0- 2 


1.67xl0- 3 


1.16xl0- 3 


9.09xl0~ 5 


3.13xl0- 5 


5.37 xlO" 5 


1.26xl0- 6 


<j{<J a ) 


(kms -1 ) 


14.27 


7.253 


4.468 


19.15 


10.48 


19.58 


16.61 


cr(ab) 


(kms -1 ) 


11.36 


5.470 


14.94 


8.877 


8.874 


26.34 


15.64 


Cr(CT c ) 


(kms -1 ) 


10.21 


2.463 


10.97 


9.450 


11.45 


29.42 


13.71 


a(a) 


(°) 


32.56 


37.80 


29.52 


25.96 


24.41 


14.58 


15.43 


a(p) 


(°) 


36.66 


35.48 


26.46 


29.98 


16.02 


33.27 


27.06 


<r(7) 


(°) 


7.58 


12.59 


22.23 


21.99 


25.78 


29.95 


17.58 



Halo local properties averaged over an ensemble of 10 spheres at each radius. The rows show 
the mass density p, the average velocity dispersion components, (<7 , <7(,, <r c ), along the principal 
axes of the velocity dispersion ellipsoid, and the angles, (a, j3, 7), between the major, intermediate, 
and minor axes of the local velocity and the global potential ellipsoids. Also listed are the 
dispersions of the above quantities. 
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Tabic 5: Semi- analytic Model Results 



Vkick 






^return 


(km s- 1 ) 


(kpc) 


(kpc) 


(Gyr) 


200 


0.0406 


0.0010 


0.0016 


280 


0.2707 


0.0010 


0.0415 


300 


0.4512 


0.0010 


0.0791 


360 


2.2022 


0.0010 


0.4735 


380 


3.7714 


0.0010 


1.6275 


400 


6.8619 


0.0010 


3.4846 


420 


10.5830 


0.0010 


8.0657 


440 


17.9090 


0.0010 


10.4097 


460 


24.0626 


0.0010 


> t H 


500 


37.2263 


1.2189 


> t H 


560 


84.6069 


1.2555 


> t H 


600 


137.3806 


17.4473 


> t u 


680 


786.7245 


276.3753 


> *H 



Columns 1,2,3, and 4 list the initial kick velocity, the MBH apocenter, its pericenter, and 
the return time within 1 pc from the center calculated using a triaxial NFW + isothermal 
spherical bulge model (see the text for details). 



